Size dependence of the minimum excitation gap in the Quantum Adiabatic Algorithm 
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We study the typical (median) value of the minimum gap in the quantum version of the Exact 
Cover problem using Quantum Monte Carlo simulations, in order to understand the complexity of 
the quantum adiabatic algorithm (QAA) for much larger sizes than before. For a range of sizes, 
N < 128, where the classical Davis-Putnam algorithm shows exponential median complexity, the 
QAA shows polynomial median complexity. The bottleneck of the algorithm is an isolated avoided 
crossing point of a Landau-Zener type (collision between the two lowest energy levels only). 
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There is considerable interest in finding optimization 
problems which could be solved much more efficiently 
by an eventual quantum computer than by a classical 
computer. An important class of classically intractable 
problems is the NP-hard category []]]. Many optimiza- 
tion problems of current interest have parameters which 
are random and so each problem corresponds to a large 
number (possibly infinite) of "instances" . The term NP- 
hard actually refers to the behavior of the computation- 
ally hardest instance, but, from a practical point of view, 
it is also of great interest to know how the time to solve 
a typical instance 0, El, the typical complexity, scales 
with problem size. Numerical studies of NP-hard prob- 
lems show that this scaling is exponential in a broad class 
of problem parameters [2J, La] • It would be a very impor- 
tant breakthrough to show that a quantum computer can 
solve the same class of problem instances of an NP-hard 
problem in less then exponential time. 

In this paper we study the typical complexity as a func- 
tion of system size for a particular quantum algorithm, 
the quantum adiabatic algorithm (QAA) proposed by 
Farhi et al. [4j . The idea is that one adds to a "problem" 
Hamiltonian, Hp, whose ground state represents a solu- 
tion of a classical optimization problem a non-commuting 
"driver" Hamiltonian, Tip,, so the total Hamiltonian is 



H(X) = (l-X)Hp, + XH P , 



(1) 



where A = A(t) is a time dependent control parameter. 
For Tip we are interested in binary optimization problems 
expressed in terms of classical Ising spins taking values 
±1, or equivalently in terms of the z-components of the 
Pauli matrices for each spin, of. The driver Hamilto- 
nian is then simply Hp, = — Y^i=i &i where of is the 
x-component Pauli matrix. 

The control parameter X(t) is at t = 0, so H—Hp>, 
which has a trivial ground state in which all 2 N basis 
states (in the o z basis) have equal amplitude. It then 
increases with t, reaching 1 at t = T (T is the runtime 
or complexity of the algorithm), at which point H=Hp. 
If the time evolution of X(t) is sufficiently slow, the pro- 



cess will be adiabatic. Hence, starting the system in the 
ground state of Hp, (all spins aligned along x), the system 
will end up in the classical ground state, which is what 
we want, with only small probability of failure. An upper 
bound for the complexity of the QAA can be given [5|, [y] , 
in terms of the eigenstates and eigenvalues of the Hamil- 
tonian, H$ n = E n <b m , 

x 2 



T>n|maxy 10 (A)|/(A£' n 



(2) 



where AE m { n corresponds to the minimum of the first 
excitation gap AE m [ n = min>, AE(X) with AE = E\ — 
E , and V n0 (X) = \dH / dX\^ n ) . Typically, matrix 
elements of H scale as a low polynomial of a number of 
spins N and the question of whether the complexity T 
depends polynomially or exponentially with N depends 
on how the minimum gap AE m i n scales with N. The 
size dependence of the minimum gap will therefore be 
the central focus of this paper. 

It is difficult to study the typical complexity of the 
QAA analytically since A*, the value of A at the min- 
imum of the gap AE(X), is different for each instance 
with fluctuations being 0{N^ 1 / 2 ), so the ensemble aver- 
aging over random instances can only be performed after 
A* has been found for each case. In the original work of 
Farhi et al. Q , the complexity of the adiabatic algorithm 
was studied numerically by direct integration in time of 
the system with Hamiltonian H. Since the size of the 
Hilbert space increases exponentially (it is of order 2 N ) 
they were limited to very small sizes, N < 20. Subse- 
quently Hogg M considered sizes up to N = 24. These 
early papers U3] found that the complexity of the al- 
gorithm scales as a roughly as A^ 2 . However, this power 
law complexity may be an artifact of the very small sizes 
studied, so it is of great interest to determine whether the 
complexity continues to be polynomial for much larger 
sizes or whether a "crossover" to exponential complexity 
is seen. To investigate this question, it is not possible 
to include all terms in the Hilbert space (as was done 
in the early work) since this becomes much too large. 
Here we use Quantum Monte Carlo (QMC) simulations, 
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with which we can study much larger sizes because only 
a sampling of the states is performed. 

There have also been QMC simulations, see e.g. Ref. [1] 
for a discussion, in which t in Eq. ([1]) is the number of 
Monte Carlo sweeps, and one estimates how the final 
excess energy (i.e. the energy above the ground state) 
varies with the total number of sweeps T. However, this 
is a "fake" dynamics, which is not necessarily represen- 
tative L 8| of the real time unitary evolution guided by 
the Schrodinger equation. Therefore the computational 
complexity of such a procedure does not necessarily cor- 
respond to that of the quantum adiabatic algorithm 

To make a comparison with the earlier work we study 
(essentially) the same model of Tip used by Farhi et 
al. 0] ■ It corresponds to an Exact Cover problem, which 
is a particular version of a Constraint Satisfaction, a com- 
monly studied problem in the NP-hard category. In Ex- 
act Cover there are N Ising spins and M "clauses" each 
of which involves three spins (chosen at random). The 
energy of a clause is zero if one spin is —1 and the other 
two are 1, otherwise the energy is 1. Thus Tip equals 

1 M 

a— 1 

+ ^ 2 ^3+^3^i+ 3 ^1^2^ 3 )i ( 3 ) 

where a± , ui and otj, are the three spins in clause a and 
the {cTi}i=o are Pauli matrices. In the absence of the 
driver Hamiltonian, the Pauli matrices can be replaced 
by classical Ising spins taking values ±1. An instance has 
a "satisfying assignment" if there is at least one choice 
for the spins where the total energy is zero. As the ra- 
tio M/N is increased, there is a phase transition where 
the number of satisfying assignments goes to zero. The 
version used by Farhi et al. considers only instances with 
a unique satisfying assignment (USA), i.e. there is only 
one state with energy 0. This has the advantage that the 
gap AE(X) is greater than zero in both limiting cases, 
Ti = Tip, and Ti = Tip, but will have a minimum at an 
intermediate value A = A*, see Fig. [T] The aim is to 
determine the size N dependence of the typical value of 
Ai^min, averaged over many instances. 

We generate instances with a USA as follows. For each 
size A, we take M clauses and prune off (i) isolated sites, 
and (ii) clauses (think of them as triangles) which are 
only connected to other clauses at one corner, since these 
give a trivial degeneracy without changing the complex- 
ity. This leaves N' sites and M' clauses. Using the stan- 
dard Davis-Putnam-Logemann-Loveland (DPLL) [9( al- 
gorithm we then see if the remaining A' sites with M' 
clauses have a USA. For each A, we choose M to max- 
imize the probability of finding a USA. Although the 
probability of finding a USA decreases exponentially with 
A, we have easily been able to find instances for A up 
to 256 and the values of M are shown in Table HI For 
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FIG. 1: (Color online) QMC results for the gap between the 
ground state and the first excited state as a function of the 
control parameter A for one instance with N = 64. The region 
around the minimum value of the gap, A_B m i n , which occurs 
at A = A*, is blown up in the inset. 



N 


16 


32 


64 


128 


192 


256 


M 


12 


23 


44 


86 


126 


166 


a 


0.7500 


0.7188 


0.6875 


0.6719 


0.6563 


0.6484 



TABLE I: For sizes N up to 256 we show values of the number 
of clauses M for which the probability of a unique satisfying 
assignment (USA), constructed as described in the text, is 
maximized. The ratio M/N is denoted by a, and is expected 
to approach the value at the quantum phase transition a c — 
0.625 [xoj] for N -> oo. For the QMC simulations we only 
used the sizes up to N — 128. 

the sizes which we will study by QMC (A < 128) the 
DPLL algorithm clearly shows exponential complexity, 
see Fig. gj 

For each instance, we use QMC to simulate the quan- 
tum system in Eqs. JT]) and ([3]) with A' spins and M' 
clauses. We simulate an effective classical model with 
Ising spins of (t) = ±1 in which r (0 < r < (3 = T _1 ) 
is imaginary time. In practice, imaginary time is dis- 
cretized into L T "time slices" each representing At = 
(3/L T of imaginary time. For, a different model, the 1-d 
Ising chain in a transverse field we have verified numeri- 
cally [lH that the scaling behavior of the energy gap [l2| 
is the same for At — > as for finite Ar, and hence it is 
plausible that a discrete Ar will work here too. 

We calculate the time-dependent correlation function 

C(r) ^EE<* + T K( r °) ) ' W 
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FIG. 2: (Color online) A log-linear plot of the median com- 
plexity of the exact cover problem using the (classical) DPLL 
algorithm as a function of N. The straight line fit works 
well demonstrating that the complexity increases exponen- 
tially with N even for quite modest sizes. This figure is for 
samples with a USA but the data for all samples (with the 
same number of clauses M) is very similar. The inset plots 
the same data on a log-log scale. The pronounced curvature 
shows that the data can not be fitted to a power law. 



with At = 1 and L T large enough that (3AE ^> 1, so the 
system is in the ground state. For r <C /3, the correlation 
function C(r) will be a sum of exponentials 



C(t) = q + J] A n exp[-(E n - E )t] 



(5) 



where the A n are constants and q, the long time limit of 
the correlation function, is determined from 



1 N ' ( 1 Lt 

1=1 \ T T =l 



(to)) 



(6) 



At large r, the sum in Eq. ([5]) is dominated by the term 
corresponding to the first excited state, (n = 1), and so 
AE can be obtained by fitting log [C(t) — q] against r for 
large r. Figure [3] shows such a fit for an instance with 
N = 128 near the minimum gap. 

We determine A_E m j n , the minimum value of the gap 
(to the first excited state), as A is varied. Fig. [1] shows 
QMC results for the gap between the ground state and 
the first excited state as a function of the control param- 
eter A for one instance with N — 64. The inset shows 
more clearly the region of the minimum gap. The gap 
is greater than zero for both A = and 1 (a property 
of this model) and is much smaller at an intermediate 
value A* in the vicinity of the quantum phase transition. 
Each instance has to be carefully monitored to find the 
minimum gap, since A* is different for each instance. 



A= 0.6375 
AE = 0.0354 
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FIG. 3: (Color online) A log-linear plot of the time dependent 
correlation function for an instance with N = 128 near the 
minimum gap. The energy gap is the negative of the slope at 
large values of r. The number of time slices was L T = 300. 
The error bars were estimated by repeating the runs many 
(typically 100) times. 



For the largest size studied, N — 128, we found that for 
some instances it was difficult to determine q accurately 
for a range of A, because the simulation was not fully 
equilibrated; the required numnber of sweeps increases 
rapidly with N. As a result, plots of C(r) — q, see Fig. [3J 
were strongly curved. In a few cases the error in the 
computed value of q was small and the problem could 
be cured by allowing q to vary slightly away from the 
computed value when doing the fits. However, we did 
not trust this procedure if the correction to q was large. 
For the remaining 13 out of 50 instances, we were able to 
provide an upper bound for the minimum gap (from the 
range of A where q was successfully computed) and this 
turned out to be less than our eventual estimate for the 
median gap. Hence we were able to obtain reliable data 
for sizes up to N = 128. However, at present we are not 
able to study much larger sizes because of the difficulty 
in determining q. 

Since we are interested in the typical minimum gap 
(among different instances), rather than the average (or 
smallest) we show in Fig. [4] the median of the minimum 
gap for N < 128. The main figure is a log-log plot, and 
the dashed line corresponds to the median AE min varying 
as AT -0 ' 73 . The pronounced curvature in the inset (log- 
linear plot) shows that the behavior is not exponential. 
The minimum gap therefore follows a power law for this 
range of sizes, implying polynomial complexity. This re- 
sult is consistent with that found by Far hi et al. Q and 



Hogg 
et al. 



for much smaller sizes (N < 20-24). Banuls 
] studied the QAA using using matrix product 



states for sizes up to A^ = 60, but their result that the 
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FIG. 4: (Color online) A log-log plot of the median of the 
minimum gap as a function of the number of bits TV up to 
TV = 128. From the satisfactory straight line fit, it is seen 
that the median Ai5 m i n decreases as a power law, TV _M with 
(j, = 0.73 ± 0.06. The number of instances is 50 except for 
TV = 64 for which it is 45. The inset shows a log-linear plot. 
The pronounced curvature shows that the behavior is not ex- 
ponential for this range of sizes, in contrast to the classical 
DPLL algorithm, data for which is shown in Fig. [2] 

complexity becomes independent of size for TV > 40 is 
surprising and quite different from ours. 

In addition to the energy gap AE(X), we also in- 



vestigated -~d 2 E /dX 2 



^Ern^\V m\ 2 /(E m - E ) 



since this gives additional information about matrix el- 
ements near the avoided crossing point A*. We deter- 
mined this from x = fo {[Hp(t)Hp(0) - (H P ) 2 ]) dr = 
— (1 — A) 2 d 2 Eq j ' d\ 2 , finding that Vo m depends on TV very 
weakly near A ~ A*. We also found that the location of 
the maximum of —d 2 E^/dX 2 coincides to a good preci- 
sion with A* , see Fig. [5j Hence the sum in the expres- 
sion for d 2 Eo/d\ 2 is dominated by its first term (m=l) 
in the vicinity of the avoided-crossing at A*, which is 
of the Landau-Zener type (collision of E\ and Eq levels 
only). This suggests that T = h\V 10 {X*)\/[e {AE min ) 2 } 
is an accurate estimate for the algorithm complexity, 
where e < 1 is an TV- independent constant. As a re- 
sult, T - TV 2 ' 1 where n = 0.73 ± 0.06. 

To conclude, by using QMC simulations we have con- 
siderably extended the range of sizes over which the com- 
plexity of the quantum adiabatic algorithm (QAA) can 
be investigated. For sizes up to TV =128, where the 
benchmark classical algorithm for satisfiability problems 
(DPLL) shows exponential complexity, the QAA shows 
polynomial behavior of the median minimum gap, and 
hence presumably polynomial behavior of the median 
complexity (contrast Fig. 2] with Fig. (2). However, our 
results for the median do not rule out the possibility 
that some instances have exponential complexity. We 
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FIG. 5: The gap AE(X) (blue), and -d 2 E /d\ 2 (red), against 
A for an instance with TV = 128. Solid lines are cubic inter- 
polations. The location of the minimum gap, A* = 0.6306, is, 
within margin of error, equal to the maximum of —d 2 Eo/dX 2 
at A = 0.6311 (both shown by vertical dashed lines). 

also found a Landau-Zener (pairwise) character of the 
avoided crossing at the minimum gap point. 
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